System biology mediated assessment of molecular mechanism for sinapic acid against breast cancer: via network pharmacology and molecular dynamic simulation

Sinapic acid is a hydroxycinnamic acid widespread in the plant kingdom, known to be a potent anti-oxidant used for the treatment of cancer, infections, oxidative stress, and inflammation. However, the mode of action for its chemotherapeutic properties has yet not been unleashed. Hence, we aimed to identify potential targets to propose a possible molecular mechanism for sinapic acid against breast cancer. We utilized multiple system biology tools and databases like DisGeNET, DIGEP-Pred, Cytoscape, STRING, AutoDock 4.2, AutoDock vina, Schrodinger, and gromacs to predict a probable molecular mechanism for sinapic acid against breast cancer. Targets for the disease breast cancer, were identified via DisGeNET database which were further matched with proteins predicted to be modulated by sinapic acid. In addition, KEGG pathway analysis was used to identify pathways; a protein-pathway network was constructed via Cytoscape. Molecular docking was performed using three different algorithms followed by molecular dynamic simulations and MMPBSA analysis. Moreover, cluster analysis and gene ontology (GO) analysis were performed. A total of 6776 targets were identified for breast cancer; 95.38% of genes predicted to be modulated by sinapic acid were common with genes of breast cancer. The ‘Pathways in cancer’ was predicted to be modulated by most umber of proteins. Further, PRKCA, CASP8, and CTNNB1 were predicted to be the top 3 hub genes. In addition, molecular docking studies revealed CYP3A4, CYP1A1, and SIRT1 to be the lead proteins identified from AutoDock 4.2, AutoDock Vina, and Schrodinger suite Glide respectively. Molecular dynamic simulation and MMPBSA were performed for the complex of sinapic acid with above mentioned proteins which revealed a stable complex throughout simulation. The predictions revealed that the mechanism of sinapic acid in breast cancer may be due to regulation of multiple proteins like CTNNB1, PRKCA, CASP8, SIRT1, and cytochrome enzymes (CYP1A1 & CYP3A4); the majorly regulated pathway was predicted to be ‘Pathways in cancer’. This indicates the rationale for sinapic acid to be used in the treatment of breast cancer. However, these are predictions and need to be validated and looked upon in-depth to confirm the exact mechanism of sinapic acid in the treatment of breast cancer; this is future scope as well as a drawback of the current study.


In-silico molecular docking
In the present study, sinapic acid was docked with 62 proteins using 3 different algorithms; AutoDock Vina, AutoDock 4.2, and Schrodinger suite Glide.Initially, the ligand was prepared by minimizing its energy, followed by target preparation and molecular docking.In addition, docking was performed on the top 5 targets with their respective standard ligands using the three algorithms.

Ligand preparation
The 3D conformation of sinapic acid was retrieved in .sdfformat from the PubChem database and was converted into .pdbformat using Discovery studio visualizer (BOVIA Discovery Studio Visualizer; https:// disco ver.3ds.com/ disco very-studio-visua lizer-downl oad) 2019.The energy of the ligand was minimized using MMFF97 forcefield 31 and was converted into .pdbqtformat for docking.

Homology modeling and target preparation
The structures of targets were initially queried in UniProt (https:// www.unipr ot.org/) database to identify available structures in Protein Data Bank (RCSB; https:// www.rcsb.org/).The targets not available were further modeled using the known FASTA sequence deposited in UniProt database using SWISS-MODEL 32 (https:// swiss model.expasy.org) (Supp.File 4 (Sheet 1)).All the hetero-atoms present in the protein were removed and saved in .pdbformat utilizing Discovery studio visualizer.

Ligand-protein docking
AutoDock vina.The ligand sinapic acid was docked against the identified proteins using AutoDock Vina at PyRx (https:// pyrx.sourc eforge.io/) ver.0.8 platform to obtain binding affinity of sinapic acid with multiple targets involved in the pathogenesis of breast cancer.Further, the targets with which sinapic acid possessed the least binding energy were visualized in Discovery Studio Visualizer and the poses of ligand pertaining to the least binding energy with maximum intermolecular interactions was selected to be foreseen in Discovery studio.AutoDock 4.2.The proteins were read in .pdbformat, previously created via Discovery studio by removing hetero atoms and unwanted chains (if any).Polar hydrogens were added along with Kollman charges to build the protein into .pdbqt.The grid box was set at the center of the macromolecule and docking was performed using a genetic algorithm with 10 genetic algorithm runs with two-point crossover mode.The pose possessing the least binding energy was visualized for protein-ligand interaction.

Schrodinger suite glide
The LigPrep module of Schrodinger suite (https:// www.schro dinger.com/) 2019 was used to prepare the ligand molecule, a low-energy conformation of the ligand was made prior to docking.Further, the protein preparation wizard module was utilized to prepare the protein by adding missing amino acids and the OPLS-3 forcefield was used to minimize the energy.Water molecules were removed and conceivable ionization states for the hetero atoms in the protein were produced; one with the highest stability was chosen.
The site map module was used to identify the binding pocket of the proteins, and a grid box was set with dimensions (18 Å, 18 Å, 18 Å) using the Glide grid generation wizard (https:// nodep it.com/ node/ com.schro dinger.knime.node.gridg en.GridG enNod eFact ory) at the center of the identified binding pocket.The Glide module of Schrodinger suite 2019 was used to perform protein-ligand docking using extra precision (XP) mode and OPLS-3 force field was implemented to dock in flexible docking mode; creating confirmations for the input ligand.Further, the Docking score was generated automatically after the completion and results were analyzed via the Glide modules XP visualizer.
Gromacs version 2022.1 (https:// www.groma cs.org/) was used to perform the MD Simulations.The topology of the protein was generated by applying charmm36 forcefield 33 using the pdb2gmx module of gromacs.The proteins were solvated using three-point water model in a dodecahedron box with dimensions of 1 nm in all directions.The model system was neutralized by adding sodium (Na + ) and Chloride (Cl -) ions as counter ions.Energy minimization was performed using steepest descent integrator with verlet cutoff scheme for a maximum of 50,000 steps to achieve the least energy confirmation.The system was equilibrated using Canonical (NVT) and Isobaric (NPT) for 100 picoseconds.V-rescale, a modified Berendsen thermostat was applied to maintain constant volume and temperature at 300 K. Similarly, a C-rescale pressure coupling algorithm was applied to maintain constant pressure at 1 bar.Particle Mesh Ewald was applied for computing long-range electrostatics, coulomb, and vander waals with a cut-off of 1.2 nm.The LINCS algorithm was used to constrain bond length.Each complex was subjected to MD run for 200 ns; the coordinates and energies were saved at every 20 picoseconds.The trajectories generated were analyzed using in-built gromacs utilities 34 .

Molecular mechanics Poisson-Boltzmann surface area (MMPBSA) analysis
The gmx_MMPBSA module 35 was used to analyze the energy contribution parameters like Vander Waals & electrostatic molecular mechanics energy, polar contribution to the salvation energy, non-polar contribution of solute-solvent interactions to the solvation energy, non-polar contribution of attractive solute-solvent interactions to the salvation energy, total gas phase molecular mechanics energy, total solvation energy, and total relative binding energy.
The MMPBSA run was performed for 82 frames from a total of ten thousand frames with an interval of 120.The Poisson Boltzmann calculations were performed using an internal PBSA solver in a sander.The MMPBSA_ ana module was used to visualize the results obtained from the gmx_MMPBSA run 36 .

Physiochemical properties and Identification of targets
The physiochemical properties of sinapic acid revealed NHBD as 2, NHBA as 5, MolLogP as 1.09, and druglikeness score of − 1.07; sinapic acid has a molecular weight of 224 g/mol and obeys Lipinski's rule of five (Fig. 1).6776 targets related to breast cancer were identified from which 45.3% of the overall targets belonged to the class of enzymes, 7.8% were identified as receptors, and 6.3% of targets were identified as a transcription factors (Fig. 2).Similarly, the proteins predicted to be modulated by sinapic acid were majorly classified as enzymes (43.5%), and receptors (8.1%) (Fig. 2).www.nature.com/scientificreports/

Gene enrichment and network analysis
The prediction of targets modulated by sinapic acid were retrieved from the DIGEP-Pred database where we identified 65 proteins possessing Pa > Pi from which 62 proteins were found to match with proteins identified for the pathogenesis of breast cancer (Fig. 3).The gene enrichment analysis predicted 95.38% of genes modulated by sinapic acid to be involved in breast cancer and the remaining 4.6% were predicted to be anti-targets (SELL, GCLM, and GSS) (Fig. 3).STRING was used to assess protein-protein interaction which was further integrated with KEGG pathway analysis where 50 pathways were identified (Supp.file 1 (sheet 1); Table 1).A network between pathways and proteins was constructed using Cytoscape based on edge count; "Pathways in cancer" was predicted to be the majorly modulated pathway.Further, PRKCA, CASP8, and CTNNB1 were predicted to be the top 3 hub genes to be regulated by sinapic acid (Fig. 4).Additionally, "cellular senescence" was identified to possess the maximum "edge betweenness" of 102 followed by "chemical carcinogenesis" with an edge betweenness of 101 (Supp.file 1 (sheet 2))."Pathways in cancer" and PRKCA displayed the highest indegree and outdegree distribution of 17 and 16 respectively.Similarly, "Human cytomegalovirus infection" was predicted to possess the highest neighborhood connectivity of 13.2 (Supp.file 1 (sheet 3)).

Gene ontology
The data for GO terms i. .The GO of the top 5 CC, MF, and BP has been represented in the form of a chord diagram (Fig. 5).The integration of GO terms with KEGG modulated proteins predicted 79% of the genes to be involved in both the GO terms as well as KEGG modulated proteins (Fig. 3).

Cluster analysis
Cluster analysis was performed via the ClueGo addon tool in Cytoscape ver.3.9.0 with CC, MF, BP, and KEGG enriched genes.The analysis displayed 1 cluster with 38 groups and 60 (96.77%) identified genes were involved in the cluster (Supp.file 3 (sheet 1)); 201 GO terms were identified with 522 connections.Additionally, 115 GO edges possessed a kappa score of 2, 163 edges possessed a kappa value of 1, and 68 edges with a kappa score above 0.8 indicating majority of the data to be reliable and accurate (Supp.file 3 (sheet 2)).Cellular response to oxidative stress (group 37) possessed the highest percent of genes/group with 36 genes ((Supp.file 3 (sheet 3)).The kappa score matrix generated by ClueGO has been depicted in (Supp.file 3 (sheet 4)).The Supp. file 5 represents cluster embraced with GO terms and KEGG enriched proteins.

Molecular stability of docked complexes
Sinapic acid-protein kinase C-α complex On completion of MD run the RMSD value was analyzed for the backbone as well as complex which was not stable at the beginning of the MD run up to 45 ns and possessed fluctuation between 1 Å to 2.8 Å.The RMSD value remained unstable till 75 ns and displayed a high peak in RMSD at 70 ns.Further, there was a drop in the RMSD after 75 ns which gradually increased up to 120 ns and became stable at 130 ns with a value of ~ 2 Å for complex and ~ 4 Å for backbone; there was a difference of ~ 2 Å observed between RMSD of the backbone and complex.The RMSF value did not display major fluctuations and was in the range of ~ 1 Å to ~ 6 Å.The radius of gyration (Rg) displayed uniformity throughout the simulation with minor fluctuations in the range of ~ 1 Å to ~ 3 Å.The number of hydrogen bonds being formed throughout the simulation was analyzed; a maximum of 4 hydrogen bonds were formed.However, there was no constant bond formation between the ligand and protein.
The solvent-accessible surface area helps to identify a free surface for the ligand to bind with protein; was found to be in the range of 35 to 45 nm 3 throughout the simulation (Fig. 8).www.nature.com/scientificreports/

Sinapic acid-caspase 8 complex
The RMSD of protein-ligand complex was found to be stabilized at ~ 105 ns where a spike increase in the RMSD was observed for backbone and complex.The RMSD fluctuation between 105 to 135 ns was ~ 3 Å thereafter it The highest number of hydrogen bonds formed between the ligand and protein was 5.However, there were no constant hydrogen bonds seen throughout the MD run.Moreover, the solvent-accessible surface area displayed a varied range of values till 100 ns simulation run; thereafter, a slight decrease in the SASA value from 40 nm 2 to 35 nm 2 was observed (Fig. 8).

Sinapic acid-β-catenin 1 complex
The RMSD of β-catenin 1-sinapic acid complex ranged from ~ 1 Å to ~ 6 Å throughout the MD run; the RMSD became stable for the backbone and complex after 100 ns of MD run thereafter the difference between the RMSD of backbone and complex was ~ 1 Å to ~ 2.5 Å.The RMSF value of the complex ranged from ~ 1 Å to ~ 5 Å, there were slight spikes visible at atoms ~ 1450 which may be due to hydrogen bonds being formed by the amino acids with sinapic acid.The radius of gyration stabilized after 50 ns and was in the range of 13 Å to 14 Å; the Rg value was observed to possess a deviation of 0.5 Å after 150 ns of MD run indicating stable protein-ligand interaction.
A maximum of 4 hydrogen bonds were formed by the protein and ligand.The number of continuous hydrogen bonds increased after 100 ns and an increase in number of hydrogen bonds was observed till 200 ns of MD run.The solvent-accessible surface area was found to be in the range of 62 nm 2 to 70 nm 2 and there were observed fluctuations in SASA ranging from 130 to 200 ns (Fig. 8).www.nature.com/scientificreports/

Sinapic cid-cytochrome 3A4 complex
The RMSD of CYP3A4-sinapic acid complex became stable after 50 ns, after which the RMSD value ranged from ~ 7 Å to ~ 10 Å throughout the simulation.There was a fluctuation in RMSD at ~ 175 ns; the difference between the RMSD of backbone and complex was observed to be ~ 1.5 Å.The RMSF value ranged between ~ 2 Å to ~ 10 Å throughout the simulation, there was increased spikes visible at atoms ~ 2100 to ~ 2400 which may be involved in the formation of bonds with the ligand.Additionally, the radius of gyration displayed fluctuation between the range ~ 19.0 Å to ~ 20.5 Å and was stable throughout the run.Moreover, a highest of 5 hydrogen bonds were formed between the protein and ligand; 3 constant hydrogen bonds were formed till 30 ns of MD run, thereafter the number of hydrogen bonds decreased to 1-2 bonds visible up to 80 ns.In addition, there was a fluctuation in the number of hydrogen bonds formed till 150 ns and no interaction was visible.Thereafter 150 ns of MD run 2-3 hydrogen bonds were visible which became constant with 3 hydrogen bonds at 190 ns.
Initially, there was a decrease in the solvent-accessible surface area till 30 ns of simulation which further increased from 140 nm 2 to 150 nm 2 .Later, the SASA value decreased as a stable complex was formed; SASA values varied throughout the simulation depending upon the stability of the complex.The SASA value ranged between 135 to 165 nm 2 throughout the simulation (Fig. 8).

Sinapic acid-cytochrome 1A1 complex
The RMSD of protein ligand complex remained stable throughout the simulation after 50 ns of stability time of the complex.The RMSD value ranged from ~ 2 Å to ~ 4.5 Å throughout the simulation with a difference of ~ 0.5 Å was observed between backbone and complex.The RMSF value ranged from ~ 1 Å to ~ 7 Å throughout the simulation; high peaks of ~ 5Å to ~ 7Å were observed at atoms ~ 1250, ~ 2000, and ~ 3000 indicating that these residues may be responsible for the hydrogen bonds formed.The radius of gyration value ranged from ~ 18.7 Å to ~ 20 Å throughout the MD run; there was a slight decrease in the Rg value after 50 ns of stabilization period of the complex thereafter the Rg value was in the range of ~ 18.8 Å to ~ 19.4 Å indicating a difference of 0.6 Å.The Rg value remained unstable after 150 ns which may be due to the breakage of bonds formed between protein and ligand.A highest of 5 hydrogen bonds were formed 2 hydrogen bonds remained constant throughout the MD run till 110 ns thereafter there only 1 hydrogen bond was constant throughout the run.The solvent-accessible surface area remained constant throughout the MD run with a range of 73 Å to 85Å.There was an increase in the SASA value after 150 ns (Fig. 8).

Sinapic acid-sirtuin 1 complex
The RMSD of sirtuin 1-sinapic acid complex was assessed to be in the range of ~ 5 Å to ~ 10 Å throughout the simulation after a stabilization period of 50 ns.There were fluctuations in the RMSD value throughout the whole MD run of 200 ns and a difference of ~ 2 Å to ~ 3 Å was observed between the backbone and complex.
The RMSF value was found to be in the range of ~ 3 Å to ~ 7 Å; a high-rise peak (~ 3 Å to ~ 7 Å) was observed between atoms ~ 1000 to ~ 1300.The radius of gyration ranged between 14Å and 17Å after the stabilization of the complex at ~ 75 ns.A highest of 4 hydrogen bonds were observed during the MD run which were not constant throughout the simulation.The solvent assessable surface area ranged between ~ 70 nm 2 to ~ 85 nm 2 ; there was an observed increase in the SASA after 140 ns of MD run which may be due to unstable hydrogen bonds formed during the MD run (Fig. 8).

Molecular mechanics Poisson-Boltzmann surface area (MMPBSA) analysis
MMPBSA analysis performed for 82 frames for all the six complexes revealed that the vander waals and electrostatic molecular mechanics energy were found to be least for β-catenin complex (− 0.11 ± 0.13 kcal/mol & − 36.49± 0.24 kcal/mol respectively) whereas, cytochrome 1A1 and sirtuin 1 complexes possessed the highest i.e. 0.58 ± 0.14 kcal/mol and − 35.23 ± 0.19 kcal/mol respectively.Moreover, total gas phase and total salvation energy were found to be lowest by β-catenin-sinapic acid and PRKCA complex i.e. 50.57± 0.61 and − 29.35 ± 0.11 kcal/mol respectively; whereas, the highest was found to be by CYP3A4 and β-catenin i.e. 54.47 ± 0.36 and − 25.74 ± 0.41 kcal/mol respectively.The total binding energy was found to be least for cytochrome 1A1 (23.92 ± 0.36 kcal/mol) indicating it to be the most stable complex (Table 2).

Discussion
The present study aimed to propose the possible molecular mechanism of sinapic acid in the treatment of breast cancer via integrating multiple system biology tools like gene set enrichment analysis, gene ontology, and molecular dynamic simulation.Sinapic acid has been reported to possess time and dose-dependent suppressive action on breast and colon cancer cell lines representing its potential to be a chemotherapeutic drug to be used for the treatment of breast cancer 37 .Also, it has been reported to possess an ameliorating effect on various organ toxicities if used in combination with chemotherapeutic drugs like cisplatin, doxorubicin, and methotrexate [24][25][26] .Sinapic acid being a plant metabolite is reported to possess ameliorating effect on various organ toxicities by chemotherapeutic agents via the regulation of Nrf2/HO-1 involved in the NF-κB signalling pathway 25 .Initially, we acquired targets involved in breast cancer through DisGeNET database where we identified 6776 genes to be involved in the pathogenesis of breast cancer, which were further matched with proteins predicted to be modulated by sinapic acid; identified from DIGEP-Pred database.These matched proteins (62 proteins) were used to construct a protein-protein interaction network via STRING.KEGG pathway analysis was used to identify various pathways being regulated by sinapic acid; 50 pathways were identified to be involved.Further, a protein-pathway network was constructed using Cytoscape ver.3.9.0 from which three proteins were identified to be lead hits (CASP8, CTNNB1, and PRKCA).GO analysis was performed based on three GO terms CC, MF, and BP.In addition, Cluster analysis was also performed to identify various groups of proteins and how they Gene ontology analysis revealed "extracellular space", "protein binding", and "response to chemical" to possess the least false discovery rate in GO analysis; CC, MF, and BP respectively.Also, we predicted that 79% of genes involved in GO and KEGG mediated proteins were in common which represents a good interaction between the GO terms and KEGG mediated proteins.Further, cluster analysis revealed the presence of "cellular response to oxidative stress" to be majorly modulated by the most number of genes (36 genes).
We performed docking on 62 targets with sinapic acid using different software's; it was predicted that cytochrome enzymes (CYP1A1, CYP3A4, and CYP3A7), CAT , SIRT1, VDR, and NOS2 possessed the lowest binding energy with the ligand sinapic acid.The molecular dynamic simulation revealed that all complexes became stable after the stabilization period (~ 50 ns to ~ 100 ns); a maximum of ~ 2.5 Å to ~ 3.0 Å difference was observed in the RMSD value of the backbone and complex.The complexes of sirtuin 1 & β-catenin with sinapic acid possessed a difference in RMSD values between backbone and complex ~ 2.5 Å which falls in an acceptable range.However, there was a fluctuation in the number of hydrogen bonds although the difference in RMSD values were within the range of ~ 3 Å, this may be due to other hydrophobic and Vander Waal interactions between the ligand and proteins 39 .Similarly, we observed in the MMPBSA that the complexes possessed low Vander Waal molecular mechanics energy indicating Vander Waal forces to be acting and making the complex stable.
Network pharmacology analysis revealed classical protein kinase C α-type (PRKCA) to be involved in majority of pathways identified in breast cancer.Protein kinase C belongs to the class of serine/threonine protein kinase family of enzymes and plays a vital role in the progression of several diseases like cancer, diabetes, autoimmune diseases, heart failure, and Parkinsonism 40 .It has been reported that the expression of PRKCA is associated with endocrine resistance and poor prognosis in ER-positive (ER +) breast tumors.In addition, expression of PRKCA is elevated in triple negative breast cancer (TNBC) patients and shown to be responsible for chemotherapy resistance and metastasis 41 .PKCα acts as an upstream regulator of FOXC2, which in turn represses the expression of p120-catenin, an important component of adherens junction that acts as the anchor for E-cadherin 42 .Hence, PRKCA could be a potent target to be suppressed for a better treatment strategy to control breast cancer.However, there is still scope to understand the role of PRKCA in breast cancer as much is still to be uncovered.Figure 9 represents the predicted probable molecular mechanism of sinapic acid for the treatment of breast cancer.
In addition, CASP8 and CTNNB1 managed to be in the top 3 lead hits identified via network pharmacology; CASP8 is found to be down-regulated in breast cancer due to promoter methylation 43 .CASP8 is an important initiator of apoptosis.Absence or down-regulation of CASP8 could cause resistance to apoptosis and is correlated with unfavorable disease outcomes, such as childhood medulloblastoma and neuroblastoma.The absence or down-regulation of CASP8 may be due to epigenetic changes 44 .Hence, promoting CASP8 expression might help well in the treatment of breast cancer.Similarly, role of CTNNB1 after birth, WNT/CTNNB1 responsive stem cells are responsible for tissue homeostasis in various organs and hyperactive WNT/CTNNB1 signaling is observed in many different human cancers 45 .The first link between WNT signaling and breast cancer was established when WNT1 was identified as a proto-oncogene capable of driving mammary tumor formation in mice 46 .However, much debate and controversy persist regarding the importance of WNT signaling for the initiation, progression, or maintenance of different breast cancer subtypes 47 .The upregulation of β-catenin and its accumulation in the nucleus promotes the transcription of genes like c-Myc & cyclin D-1 48 .Hence, CTNNB1 was predicted to be down-regulated which may be one of the mechanisms by which sinapic acid prevents the progression of cancer.
Cytochrome enzymes play a major role in the activation of drugs to activated carcinogens from which CYP1A1 enzyme plays a major role as continuous exposure to inhalation chemicals and environmental carcinogens are thought to increase the level of CYP1A1 expression in extrahepatic tissues, through the aryl hydrocarbon receptor (AhR); may be due to the tendency of CYP1A1 to metabolize carcinogens.Reports suggest that upregulation of CYP3A4 is seen in 80% of breast tumors and can be used to identify tumor response in different treatments 49 .
Moreover, catalase (CAT ) is an antioxidant enzyme that catalyzes mainly the transformation of hydrogen peroxide into water and oxygen 50 .Although catalase is frequently down-regulated in tumors the underlying mechanism remains unclear.Further, silent information regulation factor 1 (sirtuin Type 1, SIRT1), as a kind Table 2. MMPBSA analysis of hub genes & sinapic acid as a complex.All the data are presented in mean ± SEM (n = 82) and unit for each parameter is kcal/mol.∆VDWAALS Vander Waals molecular mechanics energy, ∆EEL electrostatic molecular mechanics energy, ∆EPB polar contribution to the salvation energy, ∆ENPOLAR non-polar contribution of solute-solvent interactions to the solvation energy, ∆EDISPER nonpolar contribution of attractive solute-solvent interactions to the salvation energy, ∆GGAS total gas phase molecular mechanics energy, ∆GSOLV total solvation energy, ∆GTotal total relative binding energy. of NAD + dependent class III histone deacetylation enzyme, involved in tumor proliferation, invasion, and metastasis 51 .The roles of SIRT1 in breast cancer is multifaceted depending on its substrate from upstream or downstream signaling pathway.Results have displayed that SIRT1 is significantly up-regulated in breast cancer tissues and cells, which is correlated with histological grade, tumor size, and lymph node metastasis 52 .Studies report that vitamin D has been suggested to prevent and improve the prognosis of several cancers, including breast cancer 53 ; high expression of VDR in invasive breast tumors is associated with favorable prognostic factors and a low risk of breast cancer death 54 .Hence, a high VDR expression is a positive prognostic factor which may be used to identify response of tumors in different treatments.Expression of inducible nitric oxide synthase (NOS2) has been associated with poor outcome in breast cancer 55 ; the upregulation of NOS leads to tumor angiogenesis by upregulating VEGF 56 .

Conclusion
The present study aimed to identify a possible molecular mechanism of sinapic acid (SA) in the treatment of breast cancer via the utilization of multiple system biology tools like network pharmacology, Gene ontology enrichment, and molecular dynamic simulation.The predictions revealed that the mechanism of sinapic acid in breast cancer may be due to multiple pathways and proteins like β-catenin, PRKCA, CASP8, and cytochrome enzymes (CYP1A1 and CYP3A4); the majorly regulated pathway was predicted to be "Pathways in cancer".This indicates that sinapic acid can be used in the treatment of breast cancer.However, these are predictions based on previous reports, which need to be validated and looked upon in-depth to confirm the exact molecular mechanism of sinapic acid in the treatment of breast cancer; this is future scope as well as a drawback of the current study.
Additionally, GO analysis by Pearson correlation matrix for CC; predicted the confidence interval of Pearson's r coefficient for Gene Count (GC) vs Strength & false discovery rate (FDR) to be − 0.9932 to − 0.7050 & − 0.6301 to 0.8392 with p values 0.0009 & 0.614 respectively, Strength vs GC & FDR to be − 0.9932 to − 0.7050 & − 0.7589 to 0.7471 with p values 0.0009 & 0.977 respectively, and FDR vs GC & strength to be − 0.6301 to 0.8392 & − 0.7589 to 0.7471 with p values 0.6137 & 0.9768 respectively with a sample size of 7. Similarly, for MF the confidence interval of r for GC vs Strength & FDR was predicted to be − 0.9427 to − 0.4793 & − 0.6063 to 0.4903 with p values 0.0007 & 0.787 respectively, Strength vs GC & FDR to be − 0.9427 to − 0.4793 & − 0.5649 to 0.5367 with p values 0.0007 & 0.948 respectively, and FDR vs GC & strength to be − 0.6063 to 0.4903 & − 0.5649 to 0.5367 with p values 0.7871 & 0.9476 respectively with a sample size of 13.Moreover, for BP the confidence interval of r for GC vs Strength & FDR was predicted to be − 0.8401 to − 0.7228 & − 0.4338 to − 0.1554 with p values 3.08719507014992e-036 & 0.00009 respectively, Strength vs GC & FDR to be − 0.8401 to − 0.7228 & − 0.06039 to 0.2426 with p values 3.08719507014992e-036 & 0.2334 respectively, and FDR vs GC & strength to be − 0.4338 to − 0.1554 & − 0.06039 to 0.2426 with p values 0.00008 & 0.2334 respectively with a sample size of 165 (Fig. 6).

Figure 2 .
Figure 2. (a) Classification of targets identified in the pathogenesis of breast cancer, (b) Classification of targets modulated by sinapic acid which are involved in the pathogenesis of breast cancer.

Figure 3 .
Figure 3. Venn diagram representation of (a) targets involved in breast cancer (C0678222) vs matched targets, (b) targets involved in breast cancer (C0678222) vs targets regulated by sinapic acid, (c) targets of sinapic acid vs targets involved in breast cancer (C0678222) vs anti-targets, (d) GO terms (cellular component, molecular function, and biological process) vs KEGG mediated genes.

Figure 4 .
Figure 4. Protein-pathway interaction of genes modulated by sinapic acid and pathways identified to be involved.

Figure 5 .
Figure 5. Chord diagram representation of top 5 GO terms belonging to cellular components (CC), molecular function (MF), and biological process (BP).

Figure 6 .
Figure 6.Correlation matrix analysis of GO terms (1) Cellular components, (2) Molecular function, (3) Biological process where (a) represents the correlation between Strength vs Gene count vs False discovery rate as bubble colour map, (b) represents the Pearson's r coefficient of Gene count vs Strength vs False discovery rate as a heat map.

Figure 9 .
Figure 9. Proposed possible molecular mechanism of sinapic acid in breast cancer.Indicates genes identified to be modulated by sinapic acid57 .

Table 1 .
Top 15 KEGG pathways modulated by sinapic acid against breast cancer.OGC observed gene count, BGC background gene count, FDR false discovery rate.
Where, OGC: Observed gene count; BGC: Background gene count; FDR: False discovery rateLowHigh 38long to a category biologically38.Molecular docking was performed with three different docking tools i.e.AutoDock 4.2, AutoDock Vina, and Schrodinger suite glide to attain better visibility of ligand-protein docking with matched 62 targets.